#############################################################################

#     Theis Lunge method for mediation analysis               
#     Applicable situation:
#         1.The outcome is dichotomous data;
#         2.An exposure variable and two covariates: dichotomous data;
#         2.Two parallel mediators: dichotomous data.

#############################################################################
TLmethod<-function(Data){
  Data$ATemp<-Data$A
  fitm1<- glm(M1~ATemp+C1+C2,data = Data,family = "binomial")
  fitm2<- glm(M2~ATemp+C1+C2,data = Data,family = "binomial")
  
  #Construct a new data set
  levelsOfExp <- unique(Data$A)
  Data$newID <- 1:nrow(Data)
  Data1 <- Data
  Data2 <- Data
  Data1$A1 <- levelsOfExp[1]
  Data2$A1 <- levelsOfExp[2]
  tempData <- rbind(Data1,Data2)
  
  Data1 <- tempData
  Data2 <- tempData
  Data1$A2 <- levelsOfExp[1]
  Data2$A2 <- levelsOfExp[2]
  newData <- rbind(Data1,Data2)
  
  #Compute weights
  #The weight corresponding to the first mediator
  newData$ATemp<-newData$A
  tempDir <- as.matrix(predict(fitm1,type = "response",newdata=newData))[cbind(1:nrow(newData),newData$M1)]
  
  newData$ATemp<-newData$A1
  tempIndir1 <- as.matrix(predict(fitm1,type = "response",newdata=newData))[cbind(1:nrow(newData),newData$M1)]
  newData$weight1<-tempIndir1/tempDir
  
  #The weight corresponding to the second mediator
  newData$ATemp <- newData$A
  temp <- predict(fitm2,type = "response", newdata=newData)
  tempDir <- ifelse(newData$M2, temp, 1-temp)
  
  newData$ATemp <- newData$A2
  temp <- predict(fitm2,type = "response", newdata=newData)
  tempIndir2 <- ifelse(newData$M2, temp, 1-temp)
  newData$weight2 <- tempIndir2/tempDir
  
  #The final weight
  newData$weightM <- newData$weight1 * newData$weight2
  
  #Fit a suitable model to the outcome including only exposure, auxiliary exposure variables and covariates
  require(geepack)
  newData <- newData[order(newData$newID), ]
  fitNaturalEffects <- geeglm(Y ~ factor(A) +
                                factor(A1) + factor(A2)+C1+C2,
                              data=Data, family="binomial", weights=weightM,
                              id=newData$newID)
  
  return(coef(fitNaturalEffects)[1:8])
}

